// Copyright 2012 Sergio, Raul, Dorival. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.

#ifndef GEMLAB_SHAPE_H
#define GEMLAB_SHAPE_H
#include<Eigen/StdVector>
#include<vector>

// Auxiliary
//#include "auxi/array.h"
//#include "auxi/fatal.h"
//#include "auxi/util.h"
//#include "auxi/numstreams.h"
//#include "auxi/linalg.h"
//#include "mesh/vtkcelltype.h"

namespace GemLab {

// Lin2 calculates the shape functions (S) and derivatives of shape functions (dSdR) of lin2
// elements at {r,s,t} natural coordinates. The derivatives are calculated only if derivs==true.
void Lin2(Vec_t & S, Mat_t & dSdR, double r, double s, double t, bool derivs) {
    /*
         -1     0    +1
          0-----------1-->r
    */
    S(0) = 0.5 * (1.0 - r);
    S(1) = 0.5 * (1.0 + r);

    if (!derivs) {
        return;
    }

    dSdR(0,0) = -0.5;
    dSdR(1,0) =  0.5;
}

// Lin3 calculates the shape functions (S) and derivatives of shape functions (dSdR) of lin3
// elements at {r,s,t} natural coordinates. The derivatives are calculated only if derivs==true.
void Lin3(Vec_t & S, Mat_t & dSdR, double r, double s, double t, bool derivs) {
    /*
         -1     0    +1
          0-----2-----1-->r
    */
    S(0) = 0.5 * (r * r - r);
    S(1) = 0.5 * (r * r + r);
    S(2) = 1.0 -  r * r;

    if (!derivs) {
        return;
    }

    dSdR(0,0) = r - 0.5;
    dSdR(1,0) = r + 0.5;
    dSdR(2,0) = -2.0 * r;
}

// Lin5 calculates the shape functions (S) and derivatives of shape functions (dSdR) of lin5
// elements at {r,s,t} natural coordinates. The derivatives are calculated only if derivs==true.
void Lin5(Vec_t & S, Mat_t & dSdR, double r, double s, double t, bool derivs) {
    /*
            @-----@-----@-----@-----@-> r
            0     3     2     4     1
            |           |           |
           r=-1  -1/2   r=0  1/2   r=+1
    */
    S(0) = (r-1.0)*(1.0-2.0*r)*r*(-1.0-2.0*r)/6.0;
    S(1) = (1.0-2.0*r)*r*(-1.0-2.0*r)*(1.0+r)/6.0;
    S(2) = (1.0-r)*(1.0-2.0*r)*(-1.0-2.0*r)*(-1.0-r);
    S(3) = 4.0*(1.0-r)*(1.0-2.0*r)*r*(-1.0-r)/3.0;
    S(4) = 4.0*(1.0-r)*r*(-1.0-2.0*r)*(-1.0-r)/3.0;

    if (!derivs) {
        return;
    }

    dSdR(0,0) = -((1.0-2.0*r)*(r-1.0)*r)/3.0-((-2.0*r-1.0)*(r-1.0)*r)/3.0+((-2.0*r-1.0)*(1.0-2.0*r)*r)/6.0+((-2.0*r-1.0)*(1.0-2.0*r)*(r-1.0))/6.0;
    dSdR(1,0) = -((1.0-2.0*r)*r*(r+1.0))/3.0-((-2.0*r-1.0)*r*(r+1.0))/3.0+((-2.0*r-1.0)*(1.0-2.0*r)*(r+1.0))/6.0+((-2.0*r-1.0)*(1.0-2.0*r)*r)/6.0;
    dSdR(2,0) = -2.0*(1.0-2.0*r)*(-r-1.0)*(1.0-r)-2.0*(-2.0*r-1.0)*(-r-1.0)*(1.0-r)-(-2.0*r-1.0)*(1.0-2.0*r)*(1.0-r)-(-2.0*r-1.0)*(1.0-2.0*r)*(-r-1.0);
    dSdR(3,0) = -(8.0*(-r-1.0)*(1.0-r)*r)/3.0-(4.0*(1.0-2.0*r)*(1.0-r)*r)/3.0-(4.0*(1.0-2.0*r)*(-r-1.0)*r)/3.0+(4.0*(1.0-2.0*r)*(-r-1.0)*(1.0-r))/3.0;
    dSdR(4,0) = -(8.0*(-r-1.0)*(1.0-r)*r)/3.0-(4.0*(-2.0*r-1.0)*(1.0-r)*r)/3.0-(4.0*(-2.0*r-1.0)*(-r-1.0)*r)/3.0+(4.0*(-2.0*r-1.0)*(-r-1.0)*(1.0-r))/3.0;
}

// Tri3 calculates the shape functions (S) and derivatives of shape functions (dSdR) of tri3
// elements at {r,s,t} natural coordinates. The derivatives are calculated only if derivs==true.
void Tri3(Vec_t & S, Mat_t & dSdR, double r, double s, double t, bool derivs) {
    /*      s
            |
            2, (0,1)
            | ',
            |   ',
            |     ',
            |       ',
            |         ',
            |           ',
            |             ',
            |               ',
            | (0,0)           ', (1,0)
            0-------------------1 ---- r
    */
    S(0) = 1.0 - r - s;
    S(1) = r;
    S(2) = s;

    if (!derivs) {
        return;
    }

    dSdR(0,0) = -1.0;
    dSdR(1,0) =  1.0;
    dSdR(2,0) =  0.0;

    dSdR(0,1) = -1.0;
    dSdR(1,1) =  0.0;
    dSdR(2,1) =  1.0;
}

// Tri6 calculates the shape functions (S) and derivatives of shape functions (dSdR) of tri6
// elements at {r,s,t} natural coordinates. The derivatives are calculated only if derivs==true.
void Tri6(Vec_t & S, Mat_t & dSdR, double r, double s, double t, bool derivs) {
    /*      s
            |
            2, (0,1)
            | ',
            |   ',
            |     ',
            |       ',
            5         '4
            |           ',
            |             ',
            |               ',
            | (0,0)           ', (1,0)
            0---------3---------1 ---- r
    */
    S(0) = 1.0 - (r + s) * (3.0 - 2.0 * (r + s));
    S(1) = r * (2.0 * r - 1.0);
    S(2) = s * (2.0 * s - 1.0);
    S(3) = 4.0 * r * (1.0 - (r + s));
    S(4) = 4.0 * r * s;
    S(5) = 4.0 * s * (1.0 - (r + s));

    if (!derivs) {
        return;
    }

    dSdR(0,0) = -3.0 + 4.0 * (r + s);
    dSdR(1,0) =  4.0 * r - 1.0;
    dSdR(2,0) =  0.0;
    dSdR(3,0) =  4.0 - 8.0 * r - 4.0 * s;
    dSdR(4,0) =  4.0 * s;
    dSdR(5,0) = -4.0 * s;

    dSdR(0,1) = -3.0 + 4.0*(r + s);
    dSdR(1,1) =  0.0;
    dSdR(2,1) =  4.0 * s - 1.0;
    dSdR(3,1) = -4.0 * r;
    dSdR(4,1) =  4.0 * r;
    dSdR(5,1) =  4.0 - 4.0 * r - 8.0*s;
}

// Tri15 calculates the shape functions (S) and derivatives of shape functions (dSdR) of tri15
// elements at {r,s,t} natural coordinates. The derivatives are calculated only if derivs==true.
void Tri15(Vec_t & S, Mat_t & dSdR, double r, double s, double t, bool derivs) {
    /*      s
            ^
            |
          2
            @,(0,1)
            | ',
            |   ', 9
         10 @     @,
            |  14   ',   4
          5 @    @     @
            |           ',  8
         11 @  12@   @    '@
            |       13      ',
            |(0,0)            ', (1,0)
            @----@----@----@----@  --> r
          0      6    3    7     1
    */
    double pt1 = 128.0/3.0;
    double pt2 = 32.0/3.0;
    double cc  = 1.0 - r - s;
    double t1  = r  - 0.25;
    double t2  = r  - 0.5;
    double t3  = r  - 0.75;
    double t4  = s  - 0.25;
    double t5  = s  - 0.5;
    double t6  = s  - 0.75;
    double t7  = cc - 0.25;
    double t8  = cc - 0.5;
    double t9  = cc - 0.75;
    S( 0) = pt2   * cc * t7 * t8 * t9;
    S( 1) = pt2   * r  * t1 * t2 * t3;
    S( 2) = pt2   * s  * t4 * t5 * t6;
    S( 3) = 64.0  * cc * r  * t1 * t7;
    S( 4) = 64.0  * r  * s  * t1 * t4;
    S( 5) = 64.0  * s  * cc * t4 * t7;
    S( 6) = pt1   * cc * r  * t7 * t8;
    S( 7) = pt1   * cc * r  * t1 * t2;
    S( 8) = pt1   * r  * s  * t1 * t2;
    S( 9) = pt1   * r  * s  * t4 * t5;
    S(10) = pt1   * s  * cc * t4 * t5;
    S(11) = pt1   * s  * cc * t7 * t8;
    S(12) = 128.0 * r  * s  * cc * t7;
    S(13) = 128.0 * r  * s  * t1 * cc;
    S(14) = 128.0 * r  * s  * cc * t4;

    if (!derivs) {
        return;
    }

    dSdR( 0,0) = - pt2 * (t8 * t9 * (t7 + cc) + cc * t7 * (t8 + t9));
    dSdR( 1,0) = pt2 * (t2 * t3 * (t1 + r) + r * t1 * (t3 + t2));
    dSdR( 2,0) = 0.0;
    dSdR( 3,0) = 64.0 * (cc * t7 * (t1 + r) - r * t1 * (t7 + cc));
    dSdR( 4,0) = 64.0 * s * t4 * (t1 + r);
    dSdR( 5,0) = - 64.0 * s * t4 * (t7 + cc);
    dSdR( 6,0) = pt1 * (cc * t7 * t8 - r * (t8 * (t7 + cc) + cc * t7));
    dSdR( 7,0) = pt1 * (cc * (t2 * (t1 + r) + r * t1) - r * t1 * t2);
    dSdR( 8,0) = pt1 * s * (t2 * (t1 + r) + r * t1);
    dSdR( 9,0) = pt1 * s * t4 * t5;
    dSdR(10,0) = - pt1 * s * t4 * t5;
    dSdR(11,0) = - pt1 * s * (t8 * (t7 + cc) + cc * t7);
    dSdR(12,0) = 128.0 * s * (cc * t7 - r * (t7 + cc) );
    dSdR(13,0) = 128.0 * s * (cc * (t1 + r) - r * t1);
    dSdR(14,0) = 128.0 * s * t4 * (cc - r);

    dSdR( 0,1) = - pt2 * (t8 * t9 * (t7 + cc) + cc * t7 * (t8 + t9));
    dSdR( 1,1) = 0.0;
    dSdR( 2,1) = pt2 * (t5 * t6 * (t4 + s) + s * t4 * (t6 + t5));
    dSdR( 3,1) = - 64.0 * r * t1 * (t7 + cc);
    dSdR( 4,1) = 64.0 * r * t1 * (t4 + s);
    dSdR( 5,1) = 64.0 * (cc * t7 * (t4 + s) - s * t4 * (t7 + cc));
    dSdR( 6,1) = - pt1 * r * (t8 * (t7 + cc) + cc * t7);
    dSdR( 7,1) = - pt1 * r * t1 * t2;
    dSdR( 8,1) = pt1 * r * t1 * t2;
    dSdR( 9,1) = pt1 * r * (t5 * (t4 + s) + s * t4);
    dSdR(10,1) = pt1 * ((cc * (t5 * (t4 + s) + s * t4)) - s * t4 * t5);
    dSdR(11,1) = pt1 * (cc * t7 * t8 - s * (t8 * (t7 + cc) + cc * t7));
    dSdR(12,1) = 128.0 * r * (cc * t7 - s * (cc + t7) );
    dSdR(13,1) = 128.0 * r * t1 * (cc - s);
    dSdR(14,1) = 128.0 * r * (cc * (t4 + s) - s * t4);
}

// Qua4 calculates the shape functions (S) and derivatives of shape functions (dSdR) of qua4
// elements at {r,s,t} natural coordinates. The derivatives are calculated only if derivs==true.
void Qua4(Vec_t & S, Mat_t & dSdR, double r, double s, double t, bool derivs) {
    /*
          3-----------2
          |     s     |
          |     |     |
          |     +--r  |
          |           |
          |           |
          0-----------1
    */
    S(0) = (1.0 - r - s + r * s) / 4.0;
    S(1) = (1.0 + r - s - r * s) / 4.0;
    S(2) = (1.0 + r + s + r * s) / 4.0;
    S(3) = (1.0 - r + s - r * s) / 4.0;

    if (!derivs) {
        return;
    }

    dSdR(0,0) = (-1.0 + s) / 4.0;
    dSdR(0,1) = (-1.0 + r) / 4.0;
    dSdR(1,0) = (+1.0 - s) / 4.0;
    dSdR(1,1) = (-1.0 - r) / 4.0;
    dSdR(2,0) = (+1.0 + s) / 4.0;
    dSdR(2,1) = (+1.0 + r) / 4.0;
    dSdR(3,0) = (-1.0 - s) / 4.0;
    dSdR(3,1) = (+1.0 - r) / 4.0;
}

// Qua8 calculates the shape functions (S) and derivatives of shape functions (dSdR) of qua8
// elements at {r,s,t} natural coordinates. The derivatives are calculated only if derivs==true.
void Qua8(Vec_t & S, Mat_t & dSdR, double r, double s, double t, bool derivs) {
    /*
          3-----6-----2
          |     s     |
          |     |     |
          7     +--r  5
          |           |
          |           |
          0-----4-----1
    */
    S(0) = (1.0 - r) * (1.0 - s) * (- r - s - 1.0) / 4.0;
    S(1) = (1.0 + r) * (1.0 - s) * (  r - s - 1.0) / 4.0;
    S(2) = (1.0 + r) * (1.0 + s) * (  r + s - 1.0) / 4.0;
    S(3) = (1.0 - r) * (1.0 + s) * (- r + s - 1.0) / 4.0;
    S(4) = (1.0 - s) * (1.0 - r  * r) / 2.0;
    S(5) = (1.0 + r) * (1.0 - s  * s) / 2.0;
    S(6) = (1.0 + s) * (1.0 - r  * r) / 2.0;
    S(7) = (1.0 - r) * (1.0 - s  * s) / 2.0;

    if (!derivs) {
        return;
    }

    dSdR(0,0) = - (1.0 - s) * (- r - r - s) / 4.0;
    dSdR(1,0) =   (1.0 - s) * (  r + r - s) / 4.0;
    dSdR(2,0) =   (1.0 + s) * (  r + r + s) / 4.0;
    dSdR(3,0) = - (1.0 + s) * (- r - r + s) / 4.0;
    dSdR(4,0) = - (1.0 - s) * r;
    dSdR(5,0) =   (1.0 - s  * s) / 2.0;
    dSdR(6,0) = - (1.0 + s) * r;
    dSdR(7,0) = - (1.0 - s  * s) / 2.0;

    dSdR(0,1) = - (1.0 - r) * (- s - s - r) / 4.0;
    dSdR(1,1) = - (1.0 + r) * (- s - s + r) / 4.0;
    dSdR(2,1) =   (1.0 + r) * (  s + s + r) / 4.0;
    dSdR(3,1) =   (1.0 - r) * (  s + s - r) / 4.0;
    dSdR(4,1) = - (1.0 - r  * r) / 2.0;
    dSdR(5,1) = - (1.0 + r) * s;
    dSdR(6,1) =   (1.0 - r  * r) / 2.0;
    dSdR(7,1) = - (1.0 - r) * s;
}

// Qua9 calculates the shape functions (S) and derivatives of shape functions (dSdR) of qua9
// elements at {r,s,t} natural coordinates. The derivatives are calculated only if derivs==true.
void Qua9(Vec_t & S, Mat_t & dSdR, double r, double s, double t, bool derivs) {
    /*
          3-----6-----2
          |     s     |
          |     |     |
          7     8--r  5
          |           |
          |           |
          0-----4-----1
    */
    S(0) = r * (r - 1.0) * s * (s - 1.0) / 4.0;
    S(1) = r * (r + 1.0) * s * (s - 1.0) / 4.0;
    S(2) = r * (r + 1.0) * s * (s + 1.0) / 4.0;
    S(3) = r * (r - 1.0) * s * (s + 1.0) / 4.0;

    S(4) = - (r *  r - 1.0) *  s * (s - 1.0) / 2.0;
    S(5) = -  r * (r + 1.0) * (s *  s - 1.0) / 2.0;
    S(6) = - (r *  r - 1.0) *  s * (s + 1.0) / 2.0;
    S(7) = -  r * (r - 1.0) * (s *  s - 1.0) / 2.0;

    S(8) = (r * r - 1.0) * (s * s - 1.0);

    if (!derivs) {
        return;
    }

    dSdR(0,0) = (r + r - 1.0) * s * (s - 1.0) / 4.0;
    dSdR(1,0) = (r + r + 1.0) * s * (s - 1.0) / 4.0;
    dSdR(2,0) = (r + r + 1.0) * s * (s + 1.0) / 4.0;
    dSdR(3,0) = (r + r - 1.0) * s * (s + 1.0) / 4.0;

    dSdR(0,1) = r * (r - 1.0) * (s + s - 1.0) / 4.0;
    dSdR(1,1) = r * (r + 1.0) * (s + s - 1.0) / 4.0;
    dSdR(2,1) = r * (r + 1.0) * (s + s + 1.0) / 4.0;
    dSdR(3,1) = r * (r - 1.0) * (s + s + 1.0) / 4.0;

    dSdR(4,0) = - (r + r)       *  s * (s - 1.0) / 2.0;
    dSdR(5,0) = - (r + r + 1.0) * (s *  s - 1.0) / 2.0;
    dSdR(6,0) = - (r + r)       *  s * (s + 1.0) / 2.0;
    dSdR(7,0) = - (r + r - 1.0) * (s *  s - 1.0) / 2.0;

    dSdR(4,1) = - (r *  r - 1.0) * (s + s - 1.0) / 2.0;
    dSdR(5,1) = -  r * (r + 1.0) * (s + s)       / 2.0;
    dSdR(6,1) = - (r *  r - 1.0) * (s + s + 1.0) / 2.0;
    dSdR(7,1) = -  r * (r - 1.0) * (s + s)       / 2.0;

    dSdR(8,0) = 2.0 * r * (s * s - 1.0);
    dSdR(8,1) = 2.0 * s * (r * r - 1.0);
}

// Tet4 calculates the shape functions (S) and derivatives of shape functions (dSdR) of tet4
// elements at {r,s,t} natural coordinates. The derivatives are calculated only if derivs==true.
void Tet4(Vec_t & S, Mat_t & dSdR, double r, double s, double t, bool derivs) {
    /*                    t
                          |
                          3
                         /|`.
                         ||  `,
                        / |    ',
                        | |      \
                       /  |       `.
                       |  |         `,
                      /   |           `,
                      |   |             \
                     /    |              `.
                     |    |                ',
                    /     |                  \
                    |     0.,,_               `.
                   |     /     ``'-.,,__        `.
                   |    /              ``''-.,,_  ',
                  |    /                        `` 2 ,,s
                  |  ,'                       ,.-``
                 |  ,                    _,-'`
                 ' /                 ,.'`
                | /             _.-``
                '/          ,-'`
               |/      ,.-``
               /  _,-``
              1 '`
             /
            r
    */
    S(0) = 1.0 - r - s - t;
    S(1) = r;
    S(2) = s;
    S(3) = t;

    if (!derivs) {
        return;
    }

    dSdR(0,0) = -1.0;
    dSdR(1,0) =  1.0;
    dSdR(2,0) =  0.0;
    dSdR(3,0) =  0.0;

    dSdR(0,1) = -1.0;
    dSdR(1,1) =  0.0;
    dSdR(2,1) =  1.0;
    dSdR(3,1) =  0.0;

    dSdR(0,2) = -1.0;
    dSdR(1,2) =  0.0;
    dSdR(2,2) =  0.0;
    dSdR(3,2) =  1.0;
}

// Tet10 calculates the shape functions (S) and derivatives of shape functions (dSdR) of tet10
// elements at {r,s,t} natural coordinates. The derivatives are calculated only if derivs==true.
void Tet10(Vec_t & S, Mat_t & dSdR, double r, double s, double t, bool derivs) {
    /*                    t
                          |
                          3
                         /|`.
                         ||  `,
                        / |    ',
                        | |      \
                       /  |       `.
                       |  |         `,
                      /   7            9
                      |   |             \
                     /    |              `.
                     |    |                ',
                    8     |                  \
                    |     0 ,,_               `.
                   |     /     ``'-., 6         `.
                   |    /               `''-.,,_  ',
                  |    /                        ``'2 ,,s
                  |   '                       ,.-``
                 |   4                   _,-'`
                 ' /                 ,.'`
                | /             _ 5 `
                '/          ,-'`
               |/      ,.-``
               /  _,-``
              1 '`
             /
            r
    */
    double u = 1.0 - r - s - t;
    S(0) = u * (2.0 * u - 1.0);
    S(1) = r * (2.0 * r - 1.0);
    S(2) = s * (2.0 * s - 1.0);
    S(3) = t * (2.0 * t - 1.0);
    S(4) = 4.0 * u * r;
    S(5) = 4.0 * r * s;
    S(6) = 4.0 * s * u;
    S(7) = 4.0 * u * t;
    S(8) = 4.0 * r * t;
    S(9) = 4.0 * s * t;

    if (!derivs) {
        return;
    }

    dSdR(0,0) =  4.0 * (r + s + t) - 3.0;
    dSdR(1,0) =  4.0 * r - 1.0;
    dSdR(2,0) =  0.0;
    dSdR(3,0) =  0.0;
    dSdR(4,0) =  4.0 - 8.0 * r - 4.0 * s - 4.0 * t;
    dSdR(5,0) =  4.0 * s;
    dSdR(6,0) = -4.0 * s;
    dSdR(7,0) = -4.0 * t;
    dSdR(8,0) =  4.0 * t;
    dSdR(9,0) =  0.0;

    dSdR(0,1) =  4.0 * (r + s + t) - 3.0;
    dSdR(1,1) =  0.0;
    dSdR(2,1) =  4.0 * s - 1.0;
    dSdR(3,1) =  0.0;
    dSdR(4,1) = -4.0 * r;
    dSdR(5,1) =  4.0 * r;
    dSdR(6,1) =  4.0 - 4.0 * r - 8.0 * s - 4.0 * t;
    dSdR(7,1) = -4.0 * t;
    dSdR(8,1) =  0.0;
    dSdR(9,1) =  4.0 * t;

    dSdR(0,2) =  4.0 * (r + s + t) - 3.0;
    dSdR(1,2) =  0.0;
    dSdR(2,2) =  0.0;
    dSdR(3,2) =  4.0 * t - 1.0;
    dSdR(4,2) = -4.0 * r;
    dSdR(5,2) =  0.0;
    dSdR(6,2) = -4.0 * s;
    dSdR(7,2) =  4.0 - 4.0 * r - 4.0 * s - 8.0 * t;
    dSdR(8,2) =  4.0 * r;
    dSdR(9,2) =  4.0 * s;
}

// Hex8 calculates the shape functions (S) and derivatives of shape functions (dSdR) of hex8
// elements at {r,s,t} natural coordinates. The derivatives are calculated only if derivs==true.
void Hex8(Vec_t & S, Mat_t & dSdR, double r, double s, double t, bool derivs) {
    /*
                         4________________7
                       ,'|              ,'|
                     ,'  |            ,'  |
                   ,'    |          ,'    |
                 ,'      |        ,'      |
               5'===============6'        |
               |         |      |         |
               |         |      |         |
               |         0_____ | ________3
               |       ,'       |       ,'
               |     ,'         |     ,'
               |   ,'           |   ,'
               | ,'             | ,'
               1________________2'
    */
    S(0) = (1.0-r-s+r*s-t+s*t+r*t-r*s*t) / 8.0;
    S(1) = (1.0+r-s-r*s-t+s*t-r*t+r*s*t) / 8.0;
    S(2) = (1.0+r+s+r*s-t-s*t-r*t-r*s*t) / 8.0;
    S(3) = (1.0-r+s-r*s-t-s*t+r*t+r*s*t) / 8.0;
    S(4) = (1.0-r-s+r*s+t-s*t-r*t+r*s*t) / 8.0;
    S(5) = (1.0+r-s-r*s+t-s*t+r*t-r*s*t) / 8.0;
    S(6) = (1.0+r+s+r*s+t+s*t+r*t+r*s*t) / 8.0;
    S(7) = (1.0-r+s-r*s+t+s*t-r*t-r*s*t) / 8.0;

    if (!derivs) {
        return;
    }

    dSdR(0,0) = (-1.0+s+t-s*t) / 8.0;
    dSdR(0,1) = (-1.0+r+t-r*t) / 8.0;
    dSdR(0,2) = (-1.0+r+s-r*s) / 8.0;

    dSdR(1,0) = (+1.0-s-t+s*t) / 8.0;
    dSdR(1,1) = (-1.0-r+t+r*t) / 8.0;
    dSdR(1,2) = (-1.0-r+s+r*s) / 8.0;

    dSdR(2,0) = (+1.0+s-t-s*t) / 8.0;
    dSdR(2,1) = (+1.0+r-t-r*t) / 8.0;
    dSdR(2,2) = (-1.0-r-s-r*s) / 8.0;

    dSdR(3,0) = (-1.0-s+t+s*t) / 8.0;
    dSdR(3,1) = (+1.0-r-t+r*t) / 8.0;
    dSdR(3,2) = (-1.0+r-s+r*s) / 8.0;

    dSdR(4,0) = (-1.0+s-t+s*t) / 8.0;
    dSdR(4,1) = (-1.0+r-t+r*t) / 8.0;
    dSdR(4,2) = (+1.0-r-s+r*s) / 8.0;

    dSdR(5,0) = (+1.0-s+t-s*t) / 8.0;
    dSdR(5,1) = (-1.0-r-t-r*t) / 8.0;
    dSdR(5,2) = (+1.0+r-s-r*s) / 8.0;

    dSdR(6,0) = (+1.0+s+t+s*t) / 8.0;
    dSdR(6,1) = (+1.0+r+t+r*t) / 8.0;
    dSdR(6,2) = (+1.0+r+s+r*s) / 8.0;

    dSdR(7,0) = (-1.0-s-t-s*t) / 8.0;
    dSdR(7,1) = (+1.0-r+t-r*t) / 8.0;
    dSdR(7,2) = (+1.0-r+s-r*s) / 8.0;
}

// Hex20 calculates the shape functions (S) and derivatives of shape functions (dSdR) of hex20
// elements at {r,s,t} natural coordinates. The derivatives are calculated only if derivs==true.
void Hex20(Vec_t & S, Mat_t & dSdR, double r, double s, double t, bool derivs) {
    /*
                         4_______15_______7
                       ,'|              ,'|
                    12'  |            ,'  |
                   ,'    16         ,14   |
                 ,'      |        ,'      19
               5'=====13========6'        |
               |         |      |         |
               |         |      |         |
               |         0_____ | _11_____3
              17       ,'       |       ,'
               |     8'        18     ,'
               |   ,'           |   ,10
               | ,'             | ,'
               1_______9________2'
    */
    double rp1 = 1.0+r;   double rm1 = 1.0-r;
    double sp1 = 1.0+s;   double sm1 = 1.0-s;
    double tp1 = 1.0+t;   double tm1 = 1.0-t;

    S( 0) = rm1*sm1*tm1*(-r-s-t-2) / 8.0;
    S( 1) = rp1*sm1*tm1*( r-s-t-2) / 8.0;
    S( 2) = rp1*sp1*tm1*( r+s-t-2) / 8.0;
    S( 3) = rm1*sp1*tm1*(-r+s-t-2) / 8.0;
    S( 4) = rm1*sm1*tp1*(-r-s+t-2) / 8.0;
    S( 5) = rp1*sm1*tp1*( r-s+t-2) / 8.0;
    S( 6) = rp1*sp1*tp1*( r+s+t-2) / 8.0;
    S( 7) = rm1*sp1*tp1*(-r+s+t-2) / 8.0;
    S( 8) = (1.0-r*r)*sm1*tm1 / 4.0;
    S( 9) = rp1*(1.0-s*s)*tm1 / 4.0;
    S(10) = (1.0-r*r)*sp1*tm1 / 4.0;
    S(11) = rm1*(1.0-s*s)*tm1 / 4.0;
    S(12) = (1.0-r*r)*sm1*tp1 / 4.0;
    S(13) = rp1*(1.0-s*s)*tp1 / 4.0;
    S(14) = (1.0-r*r)*sp1*tp1 / 4.0;
    S(15) = rm1*(1.0-s*s)*tp1 / 4.0;
    S(16) = rm1*sm1*(1.0-t*t) / 4.0;
    S(17) = rp1*sm1*(1.0-t*t) / 4.0;
    S(18) = rp1*sp1*(1.0-t*t) / 4.0;
    S(19) = rm1*sp1*(1.0-t*t) / 4.0;

    if (!derivs) {
        return;
    }

    dSdR( 0,0) = -0.125*sm1*tm1*(-r-s-t-2.0)-0.125*rm1*sm1*tm1;
    dSdR( 1,0) =  0.125*sm1*tm1*( r-s-t-2.0)+0.125*rp1*sm1*tm1;
    dSdR( 2,0) =  0.125*sp1*tm1*( r+s-t-2.0)+0.125*rp1*sp1*tm1;
    dSdR( 3,0) = -0.125*sp1*tm1*(-r+s-t-2.0)-0.125*rm1*sp1*tm1;
    dSdR( 4,0) = -0.125*sm1*tp1*(-r-s+t-2.0)-0.125*rm1*sm1*tp1;
    dSdR( 5,0) =  0.125*sm1*tp1*( r-s+t-2.0)+0.125*rp1*sm1*tp1;
    dSdR( 6,0) =  0.125*sp1*tp1*( r+s+t-2.0)+0.125*rp1*sp1*tp1;
    dSdR( 7,0) = -0.125*sp1*tp1*(-r+s+t-2.0)-0.125*rm1*sp1*tp1;
    dSdR( 8,0) = -0.5*r*sm1*tm1;
    dSdR( 9,0) =  0.25*(1.0-s*s)*tm1;
    dSdR(10,0) = -0.5*r*sp1*tm1;
    dSdR(11,0) = -0.25*(1.0-s*s)*tm1;
    dSdR(12,0) = -0.5*r*sm1*tp1;
    dSdR(13,0) =  0.25*(1.0-s*s)*tp1;
    dSdR(14,0) = -0.5*r*sp1  *tp1;
    dSdR(15,0) = -0.25*(1.0-s*s)*tp1;
    dSdR(16,0) = -0.25*sm1*(1.0-t*t);
    dSdR(17,0) =  0.25*sm1*(1.0-t*t);
    dSdR(18,0) =  0.25*sp1*(1.0-t*t);
    dSdR(19,0) = -0.25*sp1*(1.0-t*t);

    dSdR( 0,1) = -0.125*rm1*tm1*(-r-s-t-2.0)-0.125*rm1*sm1*tm1;
    dSdR( 1,1) = -0.125*rp1*tm1*( r-s-t-2.0)-0.125*rp1*sm1*tm1;
    dSdR( 2,1) =  0.125*rp1*tm1*( r+s-t-2.0)+0.125*rp1*sp1*tm1;
    dSdR( 3,1) =  0.125*rm1*tm1*(-r+s-t-2.0)+0.125*rm1*sp1*tm1;
    dSdR( 4,1) = -0.125*rm1*tp1*(-r-s+t-2.0)-0.125*rm1*sm1*tp1;
    dSdR( 5,1) = -0.125*rp1*tp1*( r-s+t-2.0)-0.125*rp1*sm1*tp1;
    dSdR( 6,1) =  0.125*rp1*tp1*( r+s+t-2.0)+0.125*rp1*sp1*tp1;
    dSdR( 7,1) =  0.125*rm1*tp1*(-r+s+t-2.0)+0.125*rm1*sp1*tp1;
    dSdR( 8,1) = -0.25*(1.0-r*r)*tm1;
    dSdR( 9,1) = -0.5*s*rp1*tm1;
    dSdR(10,1) =  0.25*(1.0-r*r)*tm1;
    dSdR(11,1) = -0.5*s*rm1*tm1;
    dSdR(12,1) = -0.25*(1.0-r*r)*tp1;
    dSdR(13,1) = -0.5*s*rp1*tp1;
    dSdR(14,1) =  0.25*(1.0-r*r)*tp1;
    dSdR(15,1) = -0.5*s*rm1*tp1;
    dSdR(16,1) = -0.25*rm1*(1.0-t*t);
    dSdR(17,1) = -0.25*rp1*(1.0-t*t);
    dSdR(18,1) =  0.25*rp1*(1.0-t*t);
    dSdR(19,1) =  0.25*rm1*(1.0-t*t);

    dSdR( 0,2) = -0.125*rm1*sm1*(-r-s-t-2.0)-0.125*rm1*sm1*tm1;
    dSdR( 1,2) = -0.125*rp1*sm1*( r-s-t-2.0)-0.125*rp1*sm1*tm1;
    dSdR( 2,2) = -0.125*rp1*sp1*( r+s-t-2.0)-0.125*rp1*sp1*tm1;
    dSdR( 3,2) = -0.125*rm1*sp1*(-r+s-t-2.0)-0.125*rm1*sp1*tm1;
    dSdR( 4,2) =  0.125*rm1*sm1*(-r-s+t-2.0)+0.125*rm1*sm1*tp1;
    dSdR( 5,2) =  0.125*rp1*sm1*( r-s+t-2.0)+0.125*rp1*sm1*tp1;
    dSdR( 6,2) =  0.125*rp1*sp1*( r+s+t-2.0)+0.125*rp1*sp1*tp1;
    dSdR( 7,2) =  0.125*rm1*sp1*(-r+s+t-2.0)+0.125*rm1*sp1*tp1;
    dSdR( 8,2) = -0.25*(1.0-r*r)*sm1;
    dSdR( 9,2) = -0.25*rp1*(1.0-s*s);
    dSdR(10,2) = -0.25*(1.0-r*r)*sp1;
    dSdR(11,2) = -0.25*rm1*(1.0-s*s);
    dSdR(12,2) =  0.25*(1.0-r*r)*sm1;
    dSdR(13,2) =  0.25*rp1*(1.0-s*s);
    dSdR(14,2) =  0.25*(1.0-r*r)*sp1;
    dSdR(15,2) =  0.25*rm1*(1.0-s*s);
    dSdR(16,2) = -0.5*t*rm1*sm1;
    dSdR(17,2) = -0.5*t*rp1*sm1;
    dSdR(18,2) = -0.5*t*rp1*sp1;
    dSdR(19,2) = -0.5*t*rm1*sp1;
}

typedef void (*ShpFun) (Vec_t & S, Mat_t & dSdR, double r, double s, double t, bool derivs);

int Ngeo = 14;
//                       0       1       2       3       4        5       6       7       8       9       10      11       12        13
String Geo2name[] = {"lin2", "lin3", "lin5", "tri3", "tri6", "tri15", "qua4", "qua8", "qua9", "tet4", "tet10", "hex8", "hex20", "joint"};
ShpFun Geo2func[] = { Lin2 ,  Lin3 ,  Lin5 ,  Tri3 ,  Tri6 ,  Tri15 ,  Qua4 ,  Qua8 ,  Qua9 ,  Tet4 ,  Tet10 ,  Hex8 ,  Hex20 , NULL   }; // func == shape functions
int    Geo2geob[] = {     0,      0,      1,      3,      3,       4,      6,      6,      6,      9,       9,     11,      11, -1     }; // geob == basic geometry structure, ex: tri15 => tri6
int    Geo2geof[] = {    -1,     -1,     -1,      0,      1,       2,      0,      1,      1,      3,       4,      6,       7, -1     }; // geof == face geometry code, ex: tri6 or qua8 => lin3
int    Geo2gnd [] = {     1,      1,      1,      2,      2,       2,      2,      2,      2,      3,       3,      3,       3,  0     }; // gnd == geometry ndim
int    Geo2nvc [] = {     2,      3,      5,      3,      6,      15,      4,      8,      9,      4,      10,      8,      20,  0     }; // nvc == number of vertices in cell
int    Geo2nfa [] = {     0,      0,      0,      3,      3,       3,      4,      4,      4,      4,       4,      6,       6,  0     }; // nfa == number of faces/edges in cell
int    Geo2nvf [] = {     0,      0,      0,      2,      3,       5,      2,      3,      3,      3,       6,      4,       8,  0     }; // nvf == number of vertices in face/edge of cell
int    Geo2vtk [] = {     3,     21,     -1,      5,     22,      -1,      9,     23,     -1,     10,      24,     12,      25,  2     };

enum GeoID {
    LIN2   = 0,
    LIN3   = 1,
    LIN5   = 2,
    TRI3   = 3,
    TRI6   = 4,
    TRI15  = 5,
    QUA4   = 6,
    QUA8   = 7,
    QUA9   = 8,
    TET4   = 9,
    TET10  = 10,
    HEX8   = 11,
    HEX20  = 12,
    JOINT  = 13
};

std::map<String,int> Name2geo;
int __initgeomaps() {
    Name2geo["lin2" ] = LIN2;
    Name2geo["lin3" ] = LIN3;
    Name2geo["lin5" ] = LIN5;
    Name2geo["tri3" ] = TRI3;
    Name2geo["tri6" ] = TRI6;
    Name2geo["tri15"] = TRI15;
    Name2geo["qua4" ] = QUA4;
    Name2geo["qua8" ] = QUA8;
    Name2geo["qua9" ] = QUA9;
    Name2geo["tet4" ] = TET4;
    Name2geo["tet10"] = TET10;
    Name2geo["hex8" ] = HEX8;
    Name2geo["hex20"] = HEX20;
    Name2geo["joint"] = JOINT;
}
int __initgeomaps_dummy_int = __initgeomaps();

// get geometry type
int GetGeo(int gndim, int nverts) {
    if (gndim == 1 && nverts ==  2) return LIN2;
    if (gndim == 1 && nverts ==  3) return LIN3;
    if (gndim == 1 && nverts ==  5) return LIN5;
    if (gndim == 2 && nverts ==  3) return TRI3;
    if (gndim == 2 && nverts ==  6) return TRI6;
    if (gndim == 2 && nverts == 15) return TRI15;
    if (gndim == 2 && nverts ==  4) return QUA4;
    if (gndim == 2 && nverts ==  8) return QUA8;
    if (gndim == 2 && nverts ==  9) return QUA9;
    if (gndim == 3 && nverts ==  4) return TET4;
    if (gndim == 3 && nverts == 10) return TET10;
    if (gndim == 3 && nverts ==  8) return HEX8;
    if (gndim == 3 && nverts == 20) return HEX20;
    return -1;
}

// returns natural coordinates @ nodes. [geo][{0,1,2}][nne] -- the natural coordinates of nodes
#define _ 666 // dummy value
double Geo2natc[][3/*max ndim*/][20/*max num verts*/] =
{
    { // lin2
        {-1, 1, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _ },
        { _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _ },
        { _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _ },
    },
    { // lin3
        {-1, 1, 0, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _ },
        { _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _ },
        { _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _ },
    },
    { // lin5
        {-1, 1, 0, -0.5, 0.5, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _ },
        { _, _, _,    _,   _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _ },
        { _, _, _,    _,   _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _ },
    },
    { // tri3
        {0, 1, 0, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _ },
        {0, 0, 1, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _ },
        {_, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _ },
    },
    { // tri6
        {0, 1, 0, 0.5, 0.5,   0, _, _, _, _, _, _, _, _, _, _, _, _, _ },
        {0, 0, 1,   0, 0.5, 0.5, _, _, _, _, _, _, _, _, _, _, _, _, _ },
        {_, _, _,   _,   _,   _, _, _, _, _, _, _, _, _, _, _, _, _, _ },
    },
    { // tri15
        {0, 1, 0, 0.5, 0.5,   0, 0.25, 0.75, 0.75, 0.25,    0,    0, 0.25, 0.5,  0.25, _, _, _, _, _ },
        {0, 0, 1,   0, 0.5, 0.5,    0,    0, 0.25, 0.75, 0.75, 0.25, 0.25, 0.25, 0.5 , _, _, _, _, _ },
        {_, _, _,   _,   _,   _,    _,    _,    _,    _,    _,    _,    _,    _,    _, _, _, _, _, _ },
    },
    { // qua4
        {-1,  1,  1, -1, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _ },
        {-1, -1,  1,  1, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _ },
        { _,  _,  _,  _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _ },
    },
    { // qua8
        {-1,  1,  1, -1,  0,  1,  0, -1, _, _, _, _, _, _, _, _, _, _, _, _ },
        {-1, -1,  1,  1, -1,  0,  1,  0, _, _, _, _, _, _, _, _, _, _, _, _ },
        { _,  _,  _,  _,  _,  _,  _,  _, _, _, _, _, _, _, _, _, _, _, _, _ },
    },
    { // qua9
        {-1,  1,  1, -1,  0,  1,  0, -1, 0, _, _, _, _, _, _, _, _, _, _, _ },
        {-1, -1,  1,  1, -1,  0,  1,  0, 0, _, _, _, _, _, _, _, _, _, _, _ },
        { _,  _,  _,  _,  _,  _,  _,  _, _, _, _, _, _, _, _, _, _, _, _, _ },
    },
    { // tet4
        {0, 1, 0, 0, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _ },
        {0, 0, 1, 0, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _ },
        {0, 0, 0, 1, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _ },
    },
    { // tet10
        {0,  1,  0,  0, 0.5, 0.5,   0,   0, 0.5,   0, _, _, _, _, _, _, _, _, _, _ },
        {0,  0,  1,  0,   0, 0.5, 0.5,   0,   0, 0.5, _, _, _, _, _, _, _, _, _, _ },
        {0,  0,  0,  1,   0,   0,   0, 0.5, 0.5, 0.5, _, _, _, _, _, _, _, _, _, _ },
    },
    { // hex8
        {-1,  1,  1, -1, -1,  1,  1, -1, _, _, _, _, _, _, _, _, _, _, _, _ },
        {-1, -1,  1,  1, -1, -1,  1,  1, _, _, _, _, _, _, _, _, _, _, _, _ },
        {-1, -1, -1, -1,  1,  1,  1,  1, _, _, _, _, _, _, _, _, _, _, _, _ },
    },
    { // hex20
        {-1,  1,  1, -1,  -1,  1,  1, -1,   0,  1,  0, -1,   0,  1,  0, -1,  -1,  1,  1, -1},
        {-1, -1,  1,  1,  -1, -1,  1,  1,  -1,  0,  1,  0,  -1,  0,  1,  0,  -1, -1,  1,  1},
        {-1, -1, -1, -1,   1,  1,  1,  1,  -1, -1, -1, -1,   1,  1,  1,  1,   0,  0,  0,  0},
    },
    { // joint
        { _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _ },
        { _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _ },
        { _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _ },
    },
};
#undef _ // dummy value

}; // namespace GemLab

#endif // GEMLAB_SHAPE_H
